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EXPLICIT AND IMPLICIT COMPACT HIGH-RESOLUTION 
SHOCK-CAPTURING METHODS FOR MULTIDIMENSIONAL 
EULER EQUATIONS I: FORMULATION^ 


H.C. Yee^ 

NASA Ames Research Center 
Moffett Field, CA 94035, USA 


Abstract 

Two families of explicit and implicit compact high-resolution shock-capturing methods for the 
multidimensional compressible Euler equations for fluid dynamics are constructed. Some of these 
schemes can be fourth-order accurate away from discontinuities. For the semi-discrete case their 
shock-capturing properties are of the total variation diminishing (TVD), total variation bounded 
(TVB), total variation diminishing in the mean (TVDM), essentially nonoscillatory (ENO), or 
positive type of scheme for 1-D scalar hyperbolic conservation laws and are positive schemes 
in more than one dimension. These fourth-order schemes require the same grid stencil as their 
second-order non-compact cousins. One class does not require the standard matrix inversion or a 
special numerical boundary condition treatment associated with typical compact schemes. Due to 
the construction, these schemes can be viewed as approximations to genuinely multidimensional 
schemes in the sense that they might produce less distortion in spherical type shocks and are more 
accurate in vortex type flows than schemes based purely on one-dimensional extensions. However, 
one class has a more desirable high-resolution shock-capturing property and a smaller operation 
count in 3-D than the other class. The extension of these schemes to coupled nonlinear systems 
can be accomplished using the Roe approximate Riemann solver, the generalized Steger and 
Warming flux- vector splitting or the van Leer type flux- vector splitting. Modification to existing 
high-resolution second- or third-order non-compact shock-capturing computer codes is minimal. 
High-resolution shock-capturing properties can also be achieved via a variant of the second- 
order Lax-Friedrichs numerical flux without the use of Riemann solvers for coupled nonlinear 
systems with comparable operations count to their classical shock-capturing counterparts. The 
simplest extension to viscous flows can be achieved by using the standard fourth-order compact or 
non-compact formula for the viscous terms. 


® A condensed version will appe&r in the Proceedings of the 6th International Symposium on Compu- 
tational Fluid Dynamics, Sept. 4-8, 1006, Lake Tahoe, Nevada. Pull paper was published as a NASA 
TM-110304, August 1005. Submitted to J. of Comput. Phys., August 23, 1005. 

^Senior Staff Scientist 
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I. Introduction 

Spatially high-order compact schemes have attracted much attention in recent years due to their 
narrow grid stencil and a possible enhanced accuracy over their non-compact cousins. The majority 
of these developments are aimed at wave propagation phenomena in computational aeroacoustic 
and computational electromagetics, compressible shear-layer flow and direct numerical simulation 
of turbulence. The reader is referred to Carpenter et al. (1993), Lele (1992), and Davis (1989) 
for more details. The papers by Lele and Davis discuss wave resolution and phase errors for 
linear wave propagation. Although formal extension of their schemes to nonlinear systems is 
straightforward, systematic extension of their idea to minimize phase errors and enhance wave 
resolution for coupled nonlinear systems of equations remains to be seen. Unlike the standard 
compact schemes that use symmetric compact operators, most of the recent development in 
compact methods uses asymmetric compact operators. They also require additional filtering or 
numerical dissipation for high gradient flows and generate spurious oscillations across shock 
waves and contact discontinuities even with added linear numerical dissipation. At present, 
there is no systematic extension of these asymmetric compact schemes to have high-resolution 
shock-capturing capability. Hybrid schemes using these types of methods in conjunction with 
high-resolution shock-capturing methods to enhance shock resolution were also proposed (see e.g., 
Adams & Shariff 1995). A shortcoming of this type of hybridizations is that the numerical solution 
might experience a nonsmooth transition at the switch to a different type of scheme, in addition 
to being sensitive to the choice of the numerical flux or slope limiter. For 2-D and 3-D complex 
shock wave and contact surface interactions, the switch mechanism can become less trivial. 

The motivation of the present work is to construct schemes that retain some of the unique 
properties of compact schemes and have good shock resolution without resorting to the above type 
of hybridizations. The base schemes used are compact schemes with symmetric compact operators 
for ease of extension to high-resolution shock-capturing schemes. It is anticipated that the proposed 
schemes will have a larger scope of applications than the aforementioned schemes. Independently, 
Steve Davis of Mississippi State proposed a similar idea but with different construction than the 
present work (Davis 1995). Here we define compact scheme in a broader sense than the traditional 
definition. For a desired order of accuracy, a scheme is defined as compact if its grid stencil is at 
least one grid point less than standard non-compact schemes (in each spatial direction). 

This work was prompted by the unsatisfactory resolution of second-order total variation 
diminishing (TVD) schemes used for simulating mixing layer flows with coarse and non-adaptive 
grids (see e.g., Sandham & Yee 1 989), and by the idea of Abarbanel & Kumar (1988), and Cockbum 
& Shu (1994). Abarbanel and Kumar proposed a spatially fourth-order compact scheme without 
the associated tridiagonal matrix inversion of standard compact schemes. It is computationally 
more attractive than standard compact scheme constructions. Just like schemes proposed in Lele, 
and Davis (1989), their compact scheme exhibits poor shock resolution even with added linear 
numerical dissipation. Another idea by Cockbum and Shu that the author follows is the definition 
of a local mean. It is used as a reference for introducing local limiting to avoid spurious oscillation 
while keeping the formal accuracy of a class of compact schemes. However, this so called total 
variation diminishing in the mean (TVDM) idea does not completely suppress spurious oscillations 
due to the limiting of the local mean step even for scalar hyperbolic conservation laws. 
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The objective of this work is four fold. The first objective is to modify the Abarbanel-Kumar 
compact scheme to be high-resolution at discontinuities and extend this idea to include a larger 
class of high-resolution shock-capturing schemes. A more efficient variant of the form best 
suited for multidimensional steady-state computations will be discussed. The modifications and 
improvements proposed here are expected to improve nonlinear stability and shock resolution of 
their scheme with minimum degradation at smooth regions. The second objective is to extend the 
Cockbum-Shu fourth-order TVDM scheme to include a larger class of explicit and implicit high- 
resolution schemes while maintaining the TVDM property or high-resolution at discontinuities. 
A modification of their idea and the use of different numerical fluxes are proposed to minimize 
the spurious oscillations due to the TVDM operator. The third objective is to analyze the relative 
advantages and disadvantages and the usage of these two families of schemes. The fourth objective 
is to combine the compact stencil with a variant of second-order Lax-Friedrichs numerical 
flux to increase efficiency (minimize operation counts) for combustion, thermal and chemical 
nonequilibrium flow applications. This particular form can have the option of not requiring 
Riemann solvers for coupled nonlinear systems of equations. While slightly more diffusive than 
other numerical fluxes that require Riemann solvers, the cost saving is very noticeable. One 
remedy to compensate for the slight degradation in resolution is to use a less diffusive limiter, 
and/or grid clustering and grid adaptation at high gradient and shock regions. 

Some of these proposed schemes can be fourth-order accurate away from discontinuities. For 
the semi-discrete case their shock-capturing properties are of the TVD, total variation bounded 
(TVB), TVDM, essentially nonoscillatory (ENO), or positive type (Einfeldt 1988, Liu & Lax 1995) 
for nonlinear 1-D scalar hyperbolic conservation laws and are positive schemes in more than one 
dimension. See Yee (1989) and references cited therein for background. These schemes require a 
smaller grid stencil than their non-compact cousins. Both families of fourth-order schemes require 
a grid stencil of five points for TVD and positive types of schemes and 7 points for ENO types of 
schemes in each spatial direction. On the other hand, typical grid stencils for non-compact second- 
and fourth-order high-resolution shock-capturing schemes are 5-7 and 9-11 points, respectively in 
each spatial direction. The family of schemes based on Abarbanel and Kumar does not require the 
standard matrix inversion (Ciment & Leventhal 1975, Hirsh 1975) and special numerical boundary 
treatment (Carpenter et al. 1993) associated with typical compact schemes. Due to the construction, 
both families of schemes can be viewed as approximations to genuinely multidimensional schemes 
in the sense that they might produce less distortion in spherical type shocks and are more accurate 
in vortex type flows than schemes based purely on one-dimensional extensions. The degree of 
distortion and resolution in spherical type shocks depends also on the choice of flux limiters and 
the numerical flux construction. 

For more than one dimension, the added work over the non-compact second or third-order 
counterparts only involves extra vector additions in each spatial direction. No additional 
Riemann solver or additional flux limiter computations over the second or third-order non-compact 
counterparts are involved. Most of the proposed TVDM compact forms require a 5 x 5 block 
tridiagonal matrix inversion over their non-compact counterparts. The effort to modify existing 
high-resolution non-compact shock-capturing computer codes to have the compact form for both 
families is minimal. The present modification to Cockbum and Shu’s TVDM form produces 
better shock-resolution than their Lax-Friedrichs splitting variant. One advantage of the present 
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TVDM version is that a spatially and temporally fourth-order compact form for both time-accurate 
and time-marching approaches can be readily obtained as opposed to the inability to extend the 
Abarbanel and Kumar modification in a similar manner for time-accurate computations. The 
majority of these implicit schemes are especially suited for time-marching approaches to steady- 
state numerical solutions, since higher-order spatial accuracy can be achieved with minimal effort 
and the steady states are independent of the time step. The proper choice of time discretizations 
for the proposed two families of schemes for wave propagation and computational aeroacoustic 
type of applications has not been addressed. Thus, the time discretizations discussed in this paper 
might not be optimized for the particular types of flows in the sense of wave resolution and phase 
error discussed in Tam and Webb (1993), Lele, Davis (1989), and the Workshop on Aeroacoustics 
(1995). This topic is a subject of ongoing research. 

This is the first of a two part series of papers under the same topic. This part is devoted to the 
formulation and the second part is devoted to numerical results for fluid dynamics applications. 
Section II reviews the Abarbanel-Kumar scheme and the author’s modification of their scheme. 
Section III describes the extension of their scheme to high-resolution at shocks and contact 
discontinuities. It includes the formulation with forcing or nonlinear source terms. The proof that 
one of the forms is TVD for 1-D scalar hyperbolic conservation laws is presented. Extension of 
these schemes to viscous flows is also discussed. Section IV extends the Cockbum-Shu compact 
scheme to a larger family of schemes and extends the scalar schemes to the multidimensional 
Euler equations using the various Riemann solvers. The proof that a one-parameter family of the 
explicit and implicit compact schemes is TVDM for 1-D scalar hyperbolic conservation laws will 
be included. 


n. Compact Schemes for the 3-D Euler Equations 

In vector notation, the three-dimensional (3-D) compressible time-dependent Euler Equations 
in conservation form for an equilibrium real gas can be written as 

Ut + Fx + Gy Hg = 5", (2.1a) 

where Ft = ^,F^ — ^,Gy = ^, and -ffz = ^ with the U, F, G, and H vectors given by 


' p ' 


pu 


pv 


pw 

pu 


pv} -f p 


puv 


puw 

pv 

; F^ 

puv 

; G = 

pv^ +p 

; H = 

pvw 

pw 


puw 


pvw 


pw^ +p 

- e . 


.eu +pti. 


.ev + pv ^ 


,ew + pw. 


The vector 5 = 5(®,y,z,<) can be the forcing or source term depending on the problem. The 
dependent variable U is the vector of conservative variables, and (p^u,v,w,p)'^ is the vector of 
primitive variables. Here p is the density, «, v and w are the velocity components, pu, pv, and 
pw are the *-, y- and z-components of the momentum per unit volume, p is the pressure, e = 
p[e -f («^ -1- -I- w^)/2] is the total energy per unit volume, and e is the specific internal energy. 
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For a thermally perfect gas, the equation of state is 


p — pRTy 


( 2 . 2 ) 


where R is the specific gas constant, and T is the temperature with c = e(T). For constant specific 
heats (calorically perfect gas) 


e = c„T, 


(2.3) 


where c„ is the specific heat at constant volume. 

The above flow equations are restricted to non-chemically-reacting gases. If reacting gases 
were to be included, the species continuity equations involving mass transport of chemical species 
i due to a concentration gradient in the species should be added. Thus, the scalar density function 
p becomes a vector of species mass density and the corresponding F, G, H and S are also more 
complicated leading to the increase of the vector length of U , F, G, H and 5. See Anderson 
(1989) and Park (1990) for more detail. Although the discussion is restricted to non-reacting flows, 
the form of the schemes remains the same for reacting flows. Efficient implementation of these 
schemes similar to the non-compact TVD type of schemes to reacting flows can follow the same 
procedure as in Yee & Shinn (1989) and Yee (1989). Difficulty in avoiding the wrong speed of 
propagation with discontinuous data associated with the stiff source term remains to be addressed. 
See LeVeque & Yee (1990), Yee (1989), Yee et al. (1991), Sweby & Yee (1991, 1994), Lafon & 
Yee (1991, 1992), Griffiths et al. (1992a, b), Yee & Sweby (1993), and Yee (1995) for discussion 
of this subject. Note that for equilibrium real gas and nonequilibrium flows, the form of the 
Riemann solvers and flux- vector splittings are different from the perfect gas counterparts. See Yee 
(1989) and references cited therein for these formulae. 


2.1. Background 

As discussed in the introduction, there exists many compact schemes in the literature. Here a 
compact scheme that does not require the matrix inversion associated with the standard compact 
scheme is addressed. The fourth-order in space and second-order in time Abarbanel-Kumar 
compact scheme for (2.1) with 5 = 0 takes the form 

ir”+| + 4- -t- 4- 

= - (1 - 0)[A"p,F” 4- x^VyG'' yv^H^] 

4- - X-Voyz'D^F’^ - XWo^zVyG^ - X‘Vo.yV,H^] , (2.4a) 


where A* = ^,X^ = ^,X‘ = ^ and 
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and e.g.. 


V^F=- 


VyG 


V,H 

T^Otyz 

Voyz 

'^Oxz 


T^Oxy 

Vox{Uj,k,l) 

^Oy{Uj,k,l) 

T>Oz{Uj,h,l) 


T^dx + T^dy + T^dz 

'Ddy + "^dz 

"Ddx + T^dz 
T^dx + T^dy 


Uj+i,k,i — ^Uj,h,i + Uj-i^h,h 

Uj,h+i,i - 2Uj,k,i + Uj,k-i,h 

Uj,h,i+1 - + Uj,h,i-1, 


(2.4b) 

(2.4c) 

(2.4d) 


(2.4e) 

(2.4f) 


(2.4g) 

(2.4h) 

(2.4i) 

(2.4j) 

(2.4k) 


'^dy{Fj+i,h,i) = -Fj+i.fc+i.J - 2F,+i,fc,i + (2-41) 

where Uj,h,i is the discrete approximation of U at {jAx,kAy,lAz,nAt) and ^ = |. The time 
differencing is the second-order trapezoidal formula. The extra terms that contribute to the spatieilly 
fourth-order compact scheme are the last term on the left- and right-hand side of equation (2.4a). 
Without these two terms, the implicit scheme is the classical non-compact second-order central 
difference scheme, i.e., 


i;n+i + X‘V^H^+^] 

= - (1 - ^) + X^VyG^ -1- . (2.5) 

Note that in their original paper, Abarbanel and Kumar allow 0 = 1 but their formulation is 
valid only for ^ = 1 /2. See Section 2.2 for additional discussion and for a larger family of implicit 
schemes for steady-state computations. To obtain the spatially fourth-order compact differencing, 
Abarbanel and Kumar started with (2.5) with 0 = 0 (the forward Euler time discretization) and 
three-point central difference for the convection terms. They then took a Taylor series expansion 
about (x, t/, z, <) = (jAx, kAy, /Az, nAt) and obtained a modified equation of the form 


At 

Ut -h —Utt -}-••• — — 


Ay^ 


' Ax^ „ 

Fx F XXX "h 


+ ~^^yyy + 


Az^ 

Sz jj Szzz "h 


(2.6) 



0 


or 


Z7t + -Fjb + Gy + Hz “ — 




+ ^»‘g 


yyy 




+ 


(2.7) 


To obtain a fourth order spatial differencing, they modified (2.5) with 0 = 0 by subtracting out the 
square bracket term on the right-hand side of (2.7) and used (2.1) to obtain 


t^faa Gtimm Hr 


Gyyy — 
= 




^tyy Hjtyy H^yyy 
— Utzz — -Paizz — G 

d 


’yzzt 

-F,-Gy-Hz 


(2.8a) 

(2.8b) 

(2.8c) 

(2.8d) 


The terms F*®*, Gyyy and Hzzz need only be approximated to second-order due to their 

coefficients and (see equation (2.7)). Abarbanel and Kumar approximate (2.8) at 

(j A®, X: Ay, /Az, nAf ) by 






2Ay 


2Az 


(2.9a) 


Q 

Gyyy ~ ~ ’^'^Oy{H j,h,l) ~ T^Oy 




2A« 


2Az 


(2.9b) 


Hzzz w -^^Oz(f0,fc,i) - Vaz 




2Az 


2Az 


(2.9c) 


The resulting scheme is then spatially fourth-order. Following the idea of Beam & Warming (1978) 
and assuming the homogeneous property of the Euler equations with 
one can obtain an ADI delta formulation 


= -X^(l + Voyz)V^F^ - \^(l + Vo,z)VyG" - \‘{l + Vozy)T>zH^. 




( 2 . 10 ) 


One can approximate ^(-4), ^{F), and ^(C) using the same three-point central differencing. 
Due to the delta formulation, for steady- state computations one can drop the 1 2?o a , | y and \Vftz 
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terms. See the next section for a different way of deriving the similar form with a wider family 
of implicit schemes. One can also difference ^(-4), and -^(C) by a first-order upwind 

differencing, e.g., the first-order conservative or nonconservative linearized implicit operator 
developed by the author (Yee 1986) or the first-order flux vector splitting of Steger & Warming 
(1981). 


2.2. A Modification to the Abarbanel-Kumar Implicit Scheme 
for Steady-State Computations 

The implicit scheme (2.4) or (2.10) can be used for time-accurate as well as steady-state 
computations. For time-accurate computations, the scheme (0 = 1/2) is temporally second-order. 
Observe that the terms 



+ + ®0z 




( 2 . 11 ) 


appearing in the explicit side of (2.4a) can be interpreted as added second-order numerical 
dissipation to the over all scheme. For time-accurate calculations (2.11) might have some effect 
on the smearing of shock waves depending on the procedure in solving the resulting nonlinear 
systems of algebraic difference equations. After a steady state is reached, these added second- 
order numerical dissipations vanish. This fact becomes more apparent by examining the delta 
formulation (2.10). The inherent property of (2.4a) and (2.10) carries over to the high-resolution 
modification to be discussed in Section III. The reason is that, unlike the classical way of supplying 
a linear numerical dissipation, the design principle of high-resolution shock-capturing methods 
is constructed to automatically supply the appropriate dissipation from one grid point to the 
next. Any additional terms like (2.11) would further smear the shock wave. Therefore, (2.4) 
or (2.10) supply less numerical dissipation for time-marching to steady states than time-accurate 
calculations. In addition, (2.11) contributes added computation and might degrade the diagonal 
dominant properties of the implicit operator (if first-order upwind differencing were used in (2.10) 
or other relaxation methods). To overcome the undesirable property of (2.4) for steady-state 
applications, we start with the semi-discrete form and Taylor series expand the three-point central 
spatial differencing about (jA®,&At/,/Az). Instead of (2.8), the steady part of (2.1) is used to 
approximate Gyyy and i.e., replace (2.8) with 


— 

Gyyy = 
.ffzzz = 


Gyjtm ■ffzaaj 
-f* »yy -ffzyy» 

■f’azz Gytz^ 


and (2.9) with 


G 


Gj,k+i,i - Gj^h-i^ ^ 


2Ay 


2Az 


'yyy ^ 


Oy 


Fj+ i,k,i - ^ -gj,M-n - -gj.fc.i-i 


2Ax 


2Az 


(2.12a) 

(2.12b) 

(2.12c) 


(2.13a) 

(2.13b) 
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-n 


Oz 




2Ax 


2Az 


(2.13c) 


Applying a two-parameter family of explicit and implicit time discretizations to the resulting 
semi-discrete form using (2.13) yields 


e 


+ ix;: + X^{I + 2>0,z)l>yG'”+' + V (l + Vo^y)V,H^+^ 


-u- 


X^{l+Voyz)V^F^ + X^I + X>o,.)Py(?" -I- A*(J -f Vo.y)V,H^ 


+ 


0 ) 


1 -f a> 


nk,i - 


(2.14a) 


Here 0 < ^ < 1. The scheme is temporally second-order if 0 = a; -f- 1 and first-order otherwise. 
When 0 / 0, the method is implicit. When a; = 0, it recovers the one-parameter family case. That 
is, this two-parameter family includes the first-order implicit Euler (^ = 1, w = 0) and the second- 
order three-point backward differentiation (0 = 1, a> = 1/2) methods. Higher-order implicit 
discretizations can be achieved using a three-parameter family of linear multistep methods. See 
Lambert (1973) for the formula. Various iterative, preconitioning and/or relaxation methods (Saad 
1994, Turkel 1993) can be used to solve (2.14a) for time-accurate or steady-state computations. 
Note that unless 5 0 in (2.1), one cannot achieve the compact property for (2.14a) in 1-D 

because the 1-D form collapses to the standard second-order case. 

The analog of (2. 10) for the delta formulation of (2. 14a) can be readily obtained. For steady- state 
applications, the terms |(l?oy X>oz), |(7?o» + X>oz) and |(I>o» + T^oy) on the implicit left-hand 
side can be dropped. For the delta formulation, it yields 


(2,14b) 

Again, for steady-state computations, one can difference ^(A), ^(5), and ^(C') by first-order 
upwind differencing as discussed before. 

Although scheme (2.14a) can be used for time-accurate computations, the spatial accuracy is 
no longer fourth-order. In this case, it appears that (2.14a) has an added advantage over (2.5). The 
extra cross derivative terms |(2?oy + 2?oz)7^»-f’”, i(^o» + l>oz)2^yG" and ^(Vom + Voy)VzH^ 
can be viewed as approximations to genuinely multidimensional schemes in the sense that these 
terms would produce less distortion in spherical type shocks and are more accurate in vortex type 


r 

1 




1+ (I) 


+ 


\+u) dx J[ dy J[ 1 + u) dz 

A*(/-moyz)2?»i^” + A*'(l-h 2?o»z)2>y(?” + A'(l-f 2?o.y)75z.ff” 
Ki,,i - Kua 
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flows than schemes based purely on one-dimenslonal extensions. Of course the degree of distortion 
also depends on the grid, the form of the numerical flux and the flux limiters to be discussed in a 
later section. 

Using the one-leg formulation of Dahlquist (1979), an alternative to (2.14a,b) is 






P'm,. - 


(2.14c) 


with TJ = - 1 - For 0 = 1/2 and a> = 0, this one-leg formula is the well known 

mid-point implicit method. Note that the noniterative linearized form (Beam & Warming 1978, 
Yee & Sweby 1994) of the midpoint implicit formula reduces to the regular noniterative linearized 
trapezoidal formula. Also, there is no one-leg version of the backward Euler method since (2.14c) 
reduces to (2. 14a) for 0 = land a; = 0. Higherthan second-order implicit counterparts of (2. 14a,b) 
and the one-leg formulations can also be obtained in a similar manner but the resulting scheme 
involves more than three time levels. 

Iterative and/or relaxation procedures can be used to solve (2.14a) and (2.14c). If iterative 
relaxation procedures are used to solve (2.14), (2.14c) requires fewer flux evaluations and flux 
additions than (2.14a) for 0 0 and 0 7 ^ 1, and a> / 0. In this case, the linearized Jacobian of the 

fluxes at the “n” time-level for Newton-type iterative procedures can be used. 


2.3. Compact Explicit Schemes for Steady-State Computations 

One of the easiest procedures for obtaining higher than second-order time discretizations for 
multidimensional problems is the Runge-Kutta method. There are many variants of the Runge- 
Kutta method in the literature. See Lambert (1973), Butcher (1987) and Carpenter & Kennedy 
(1994) for details. The standard fourth-order Runge-Kutta takes the form 


jfei 

^2 






iZ°(t/”) 


R° (v" + 

(u'^ -F Atki^ 


R° 

R° 


Trn At 


- 1 - Atk$ 

k% -{• 2^2 ”1“ 2A/3 -|- k^ 


(2.15a) 


Shu’s third-order Runge-Kutta form that is compatible with TVD, TVB and ENO schemes takes 
the form 
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uw = u^- AtR°{U^) 

UW = \ jjr^ + + lAtR°(U^^h (2.15b) 

4 4 4 

= lu^ + ^ir(2) + ^AtR°{U^^^). 

o o u 

The four stage Runge-Kutta method proposed by Jameson et al. (1981) (fourth-order for constant 
coefficients and second-order for nonlinear problems) takes the form 

C^(i) = CT” - \AtR°(U^) 

4 

cr(2) = u^~ lAfi2°(lf(i)) 

c^(’) = u^~ -AtR°{U^^^) (2.15c) 

2 

ur^+i = t;-n _ Afi2°(tr(’>). 

Additional explicit and implicit Runge-Kutta-type methods can be found in Butcher, and Carpenter 
& Kennedy (1994). The proper choice of time discretization that is compatible with a chosen 
spatial discretizations is crucial in achieving low phase and amplitude errors for time-accurate 
computations. This subject is ongoing research. Following the same form as (2.14a), a spatially 
fourth-order compact scheme can be obtained by defining R° as (dividing the square bracket of 
right-hand-side of (2.14a) by At and setting 0 = 0 and « = 0). 


R° 


1 

A^ 


" 

1 

m m 


1 

“ • 

I + T^Oyz 

- — 

Ay 

I + ®0»z 

VyG^- 

Az 

I + 2^0»y 


(2.15d) 


Again, one can use (2.15) for time-accurate computations. The comments discussed in the 
paragraph above (2.14c) of Section 2.2 hold true for (2.15). 


m. Compact High-Resolution Shock-Capturing Schemes 

In the following discussion, wherever there is no confusion, the terms TVD, TVB, ENO or 
positive scheme are loosely used for schemes that are TVD, TVB, ENO or positive for (a) the fully 
discretized form, (b) the semidiscretized form, or (c) the frozen constant coefficient case. Note 
that all TVD, TVB and ENO schemes are a subclass of positive schemes and all TVD schemes 
are a subclass of TVB and ENO schemes. Also when we use the terms TVD, TVB or ENO, they 
mean the form has these types of properties for 1-D scalar constant coefficient hyperbolic PDEs 
or nonlinear scalar conservative laws. It is remarked that regardless of the type of high-resolution 
method for 1-D, the final scheme for multidimensions is only of the positive type of scheme in 
the sense of Einfeldt (1988), and Liu and Lax (1995). Strictly speaking, higher than first-order 
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TVD-type schemes exist only for 1-D scalar hyperbolic conservation laws and for 1-D linear 
hyperbolic systems. 

A careful examination of (2.4), (2.10), (2.14) or (2.15) reveals that spurious oscillations across 
discontinuities can be avoided if one replaces all of the three-point central differences of the fluxes 
(2.4b,c,d) by one of the spatially non-compact second-order high-resolution shock-capturing TVD, 
TVB, ENO or positive-type schemes. Note that spatially higher than second-order high-resolution 
non-compact schemes can also be used but additional analysis is needed on the over all order of 
accuracy of the final scheme. In other words, redefine (2.4b,c,d) by 


V^F = 




(3.1a) 


VyG = 




U 


(3.1b) 


V^H = 




(3.1c) 


where and are the non-compact second-order numerical fluxes to be 

defined shortly. Most of these numerical fluxes can be viewed as a spatially three-point central 
differencing with a nonlinear numerical dissipation as described in Harten (1983), Yee & Harten 
(1985) and Yee (1985b). The majority of these numerical fluxes have a 5-point grid stencil in 
each spatial direction. The attractive property of these compact high-resolution shock-capturing 
schemes is that fourth-order accuracy can be achieved using the same grid stencil and numerical 
fluxes as their second-order non-compact cousins. Using (3.1), (2.4a) becomes 




AS' 

(i-e)i + Vo^, 

/nrn 

A' 

(1 — 6)1 + 

frn jTn 




(3.2a) 


In symbolic notation 


= iJi . CT”. 


(3.2b) 


For steady-state computations (see Section 2.2), one can drop the last term on the implicit and 
explicit sides of (3.2a). A two-parameter family analog of (2.14a) is 
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+ v(i + p..,)[ir;+j^.-j;+;_.]} 

+ A‘ (I + Do.,) j ] I + irp;; 


or in symbolic notation 



(3.3b) 

Similarly one can obtain the corresponding high-resolution shock-capturing form of (2.14b), 
(2.14c) and (2.15) (i.e., with the appropriate numerical fluxes (3.1) instead of the three-point 
centred differences of (2.4b) - (2.4d)). Denote these high-resolution analogs of (2.14b), (2.14c) and 
(2.15) in symbolic notation respectively as 

Ld ■ AEfj;+; = ■ P" + 1 " „ [o-AV . 

(3.4) 


(3.5) 


(3.6) 


Here, Ld and Rd are the high-resolution analogs of the implicit and explicit operators for (2.14b), 
Ro is the analog of the one-leg operator for (2.14c) and Re is the analog of the symbolic notation 
for the multistage Runge-Kutta method for (2.5). The vector U in (3.5) is the same as in (2.14c). 

Again, for 0 7 ^ 0 and 0^1, and w 7 ^ 0 using iterative relaxation methods to solve the one-leg 
high-resolution formulation (3.5) requires fewer numerical flux evaluations and numerical flux 
additions than (3.3) and (3.4). Extensive numerical experimentation is needed to determined the 
relative convergence rate of (3.3), (3.4) and (3.5) for time-marching to the steady state. 

Variants of the established form of the second- and third-order numerical fluxes Fj+^,k,i that 
exhibit high-resolution shock-capturing capability have flooded the literature in the past six years. 
Most of these later numerical fluxes are applicable for the proposed scheme. Here, only a few of 
the established forms of the numerical fluxes are given. 

For the Harten and Yee upwind TVD-type scheme and Yee-Roe-Davis symmetric TVD-type 
schemes, the numerical flux using the Roe’s approximate Riemann solver is of the form (Yee & 
Harten 1985, Yee 1985b) 
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Pj+i,k,i + Fj,k,i + Rj+ 



(3.7) 


Here -Rj+ » is the right eigenvector matrix of using Roe’s approximate average state and 
has several forms which are now discussed separately. 


Second- Order Symmetric TVD Scheme . The elements of the vector ^ 1 , denoted by ( i ) , 
for a general second-order symmetric TVD scheme are 

- <?y+5l- (3-8») 


The value a* ^ is the characteristic speed where I = 1,2,..., 5, of evaluated at some 
symmetric average of Uj^h.i and Z7j+i,fc,z. The function ^ is an entropy correction to 1 1. One 
possible form is (Harten & Hyman 1983) 





[( 


a . 1 

i+i 


f + 




(3.8b) 


a* is the Ith jump in the characteristic variable in the aj-direction. For problems containing only 
unsteady shocks. Si is usually set to zero. Note that entropy-violating phenomena occur only for 
steady or nearly steady shocks. For steady-state problems containing strong shock waves, a proper 
control of the size of is very important, especially for hypersonic blunt-body flows. See Yee 
et al. (1991) for a discussion. The ‘limiter’ function Q\, i , expressed in terms of the Jump in the 
characteristic variables, can be of the form 


Qj+t = minmod(a'._, ) -f minmod(a*.^» ,aj.^») - (3.8c) 

= minmod(a*._i,a*.^i,o5.^»), (3.8d) 

Qj+i =minmod[2a'_2,2a'^a,2a5.^a,i(a'_x -haj.+ s)]. (3.8e) 


The minmod function of a list of arguments is equal to the smallest number in absolute value if the 
list of arguments is of the same sign, or is equal to zero if any arguments are of opposite sign. 


Second- Order Upwind TVD Scheme : The elements of the vector denoted by (^‘-^i)^ 
for a second-order upwind TVD scheme, originally developed by Harten and later modified and 
generalized by Yee (1985b), are 



(3.9a) 
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Examples of the ‘limiter’ function gj can be expressed as 


(3.9b) 


g] = minmod(aJ._i,aJ.j.i), 


(3.9c) 

(3.9d) 


'] = 1 [(«*•+! f + ^2] + ot 5 + 1 )* + ^2] }/ [("jf+i) + ^ 

I 

g] = minmod(2a5._,,2aJ.^i, ^(a' + i + 

g) = S ■ max 0, min(2|a*.^ 1 1, 5 • a'._ 1 ), min(|aj.^i |, 25 • aJ-_ 1 ) ; 5 = sgn(a'.^ 1 ). 


9e) 
(3.9f) 


(3.9g) 


Here 82 is a small parameter to prevent division by zero and sgn(a*.^ 1 ) = sign(a*.^ 1 )• In practical 
calculations lO"^ < 82 < 10“® is a commonly used range. For o*.,.! + ~ d] is set to 

zero in (3.9d). See Sweby (1984) for the construction of additional limiter functions. 

Later development in limiters has flooded the literature and has created much debate. Most 
of the improvements are usually problem dependent. See Donat (1994) on the error propagation 
for nonlinear approximations to hyperbolic equations containing discontinuities in derivatives. 
For the last six years of development in flux limiters, see articles which have appeared in the 
Journal of Computational Physics, the International Journal for Numerical Methods in Fluids, the 
Journal of Computers & Fluids, the Proceedings of the AIAA conference in CFD, the Proceedings 
of the International conference in Numerical Methods for Fluid D 5 mamics, the Proceedings of 
the Symposium for Computational Fluid Dynamics, the Proceedings of the GAMM-Conference 
on Numerical Methods in Fluid Mechanics, and the Proceedings of the conferences organized 
by The Institute for Computational Fluid Dynamics of the Universities of Oxford and Reading, 
England. For high gradient and/or high frequency wave propagation with shock waves and 
aeroacoustics applications, suitable limiters and the proper amount of limiting are essential to the 
overall accuracy of flow computations in addition to the comments above Section III. See Davis 
(1995) for a possible limiter for this type of flow. Extensive numerical experimentation is needed 
to determine the performance of these schemes and limiter combinations, and will be reported in a 
forthcoming paper. 


Shu (1987, 1988) showed procedures for modifying some existing TVD schemes such that the 
resulting schemes could be proven to be TVB and of globally higher-order accuracy in space, 
including extrema points. For example, by replacing the g\ function by gf^ as discussed in Shu, 
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the modified flux schemes can be made uniformly second-order accurate even at points of extrema. 
One of the forms suggested by Shu is 

gf =^minmod(aJ.^i,u>aj._| -1- JWA®2sgn(a*.^,)) 

-|-^minmod(a*._i,wa‘.^i -1- Af A®^sgn(a*._ , )). (3.10) 

Herel < o> < 3andM > 0. Shu suggests setting M = 50 for the Burgers’ equation computations. 

One can also modify the Osher-Chakravarthy (1986) method by changing the appropriate flux 
limiter functions as above. Modification of other methods can be found in Shu (1988, 1989). 


High-Resolution TVD, TVS & ENO Lax- Friedrichs Schemes: The corresponding high- 

resolution TVD Lax-Friedrichs schemes for system cases can be obtained by defining the ip 

^mam 

function to be (see Shu (1987, 1988) ) for any of the {<f>K i )^ or (^*. , i )^. The 

value can be | | + ) where u is the velocity in the «-direction and c is the sound 

speed. See Shu (1989) for additional formulae. In addition by redefining the ip, one can obtain a 
high-resolution TVB Lax-Friedrichs method by changing the limiter function to be the appropriate 
form as discussed above. Although using the Lax-Friedrichs numerical flux would introduce more 
numerical dissipation into the scheme, an entropy inequality is automatically satisfied with this 
numerical flux. Thus one does not have to deal with an arbitrary parameter ^i. In addition, at each 
grid point a savings of a (5 x 5) matrix- vector multiplication is realized in each direction. 

High-Resolution ENO and Positive Schemes : If ENO (Harten & Osher 1987) and positive 
schemes other than the Lax-Friedrichs numerical flux discussed in the previous paragraph are 
desired, one can define the appropriate numerical flux according to the references cited earlier. 


3.1. The MUSCL Approach 


MUSCL Approach Using an Approximate Riemann Solver: The numerical flux function 

for an upwind MUSCL-type scheme as described in Yee (1985a, 1989) using the 
local-characteristic approach can be expressed as 


The elements of and the vector (o£°)j^. i are given by 


(3.11a) 



(3.11b) 

(3.11c) 
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where ^{{a° i ) can be |(a° 1 1 or the same form as (3.8b). Here (a° i are the eigenvalues 

and is the eigenvector matrix of evaluated using a symmetric average between 
and U^,t\ i.e., 

3+S 


(«»)'+. =a'(tf* (3.11d) 

Rj+i (3.1U) 

However, there are options in applying the limiters for system cases. Namely, one can impose 
limiters on the conservative, primitive, or characteristic variables. 

Various “slope” limiters can be used to eliminate unwanted oscillations. A popular one is the 
“minmod” limiter; it modifies the upwind-biased interpolation as follows: 

j [(1 - + (1 + 5 ) A ^] , 

pA . = j [(1 - v)^, + (1 + 5)A^1 , 

= minmod(A,.+ i,a;A,._i), 

= minmod(Aj-+i,a;Aj-+|), 

where 

minmod(p, a>g) = sgn(p) • max|o,min[|pl,a>gsgn(p)] j 

with p = Aj^i and q = Aj_ i in equation (3.12c). Here the spatial order of accuracy (before the 
application of limiters) is determined by the value of ij: 

Tj = —1, fully upwind scheme 
^ = 0, Fromm scheme 
^ = 1/3, third-order upwind-biased scheme 
^ = 1, three-point central-difference scheme 

and 1 < o> < with 

One can improve the global order of accuracy (TVB) of the MUSCL scheme (3. 1 1) by modifying 
17^1 and in equations (3.12) by 


(3.12a) 

(3.12b) 

(3.12c) 

(3.12d) 

(3.12e) 

(3.12f) 
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where 


- 

= Uj+i,h,i - 4 [(1 - v)^j+* + (1 + v)^j+^ ] y 



^Af 

Aj.^i = minmod(Aj.^i,a>Aj._i + JWA®^sgn(Aj-^i )), 

- — w-M 

Aj-+i = minmod(Aj-^i,a;Aj-^| + M A*^sgn(Aj-^i )), 


(3.13a) 

(3.13b) 

(3.13c) 

(3.13d) 


minmod(p, 5) = sgn(p) • tnax| 0, min [|p|, gsgn(p)] | , (3.13e) 

withp = Aj-^i and q = wAj_i + AfAa^sgn(Aj_,.i ) in equation (3.13c). 

MUSCL Approach Using the Lax- Friedrichs Numerical Fluxr. The numerical flux function 
j for a MUSCL-type approach using the Lax-Friedrichs numerical flux can be expressed as 

^,+},M= S W+i) + nP-?n) + *,-+i]. (3.14.) 


where , 1 is 
j+j 


and 


*5+1 


( o \maa 

(JJ^ TT^ \ 


(3.14b) 


K)r;i = o(Kv^i+sVi)- 


(3.14c) 


There is a tremendous savings in operation count (especially for multidimensional problems 
and/or nonequilibrium flows) in using the MUSQ^-Lax-Friedrichs numerical flux instead of the 
Roe-type first-order upwind numerical flux when the limiter function is applied to the conservative 
or primitive variables instead of the characteristic variables. In problems containing contact 
discontinuities as well as shocks, one can use a more compressive limiter for the density and a less 
compressive limiter for the other variables. Note that one does not have a similar savings using the 
Lax-Friedrichs numerical flux for the non-MUSCL formulations. 


MUSCL Approach Using Flux-Vector SpliUings : The numerical flux 1 ,fc,j for either flux- 

vector splittings, can be expressed as 




(3.15) 
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where ) are evaluated using either the generalized Steger-Warming splitting or the 

generalized van Leer splitting. The vectors and are the same as in equation (3.12). See 
Yee (1989) and references cited therein. 

3.2. Compact High-Resolution Schemes for Problems Containing Source Terms 

When S in (2.1) is not zero and 5 is a function of U, x,y and z, if one uses a pointwise 
evaluation of S{U,x,y,z), i.e., S{U,x,y,z) « S{Uj,h,i,j^x,kAyJAz), (2.7) becomes 


Ut + F^ + Gy + H,-S = 


'At„ Ax^ „ 
2 3! 




yyv 


Az^ 


+ ••• , (3.16) 


and (2.8) becomes 


i^a*a = iS'aa — 17fa* — “ -ffzaa (3.17a) 

^yyy ~ ^yy ^^yy -®^*yy (3.17b) 

Hzzz — Szz U^tzz F^zz Gyzz (3.17c) 

= ^\s-F^-Gy-H}. (3.17d) 


For the semi-discrete formulation for steady- state computations, one drops all the terms containing 
the t-derivatives. Thus the analogy for (3.2), (3.3) and (3.5) for problems containing source terms, 
respectively, are 


i. ■ £?’•+> - 9A(SJ+} = fl. ■ jr" + A<|(1 - e)I + , (3.18) 




UL^t . t U) 

1 -t- a> 1 -h o> 


7T” rrn-1 


+ - — I — A<2?( 
1 + a; 


5 n 

(3-19) 


nt', = ^ + 1^ - t'w 


1-9 

+ j-^p-^At2?0ayz<5j,fc,/ 


(3.20) 


The analogy for (3.4) and (3.6) can be obtained in a similar manner. 


From the studies in LeVeque & Yee (1990), Griffiths et al. (1992), and Lafon & Yee (1991, 
1992), pointwise evaluation might not be the optimal discretization for the source terms (in terms 
of stability and accuracy), depending on the method of discretizing the convection term. Readers 
are referred to these references for additional discussion. 
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3.3. TVD Property for the scalar 1-D Case 

For 1-D, all the compact form collapses to the standard second-order non-compact case unless 
5 ^ 0 in (2.1a). For (3.2), the 1-D version has two extra terms and |2?o*(«^) over 

the standard second-order non-compact case. Therefore, only the TVD property for the scalar 1-D 
form of (3.2) will be discussed. Consider a 1-D hyperbolic conservation law 

ut -f /a = 0. (3.21) 

A one-parameter family of five-point two time-level 1-D version of (3.2) for (3.21) with A = ^ 
is given by 


u 


n+1 


+ - /;+|) + - ^(1 - i i (3.22) 


fn+l^ 


,"+l^ _ 


where 0 < 0 < 1, the numerical flux f?^i = /(«"_i>«”>«”+u«"+ 2 )> fj^i = 
/(«"![■/, For the proof later, rewrite X>oa(«i) as 

2>o»(«j) = («i+i - «i) - («i - uj-i) = 1 - A_,._ 1 . (3.23) 

To simplify the notation, rewrite equation (3.22) as 


L . = i? . u". 

The total variation of a mesh function u" is defined to be 


Here the shorthand notation 


OO 

TF(«")= X) 

j=-oo 


(3.24) 


(3.25) 


Aj+ 1 = «j+x - «i (3.26) 

for any mesh function u is used. The numerical scheme (3.22), for the initial- value problem (3.21) 
is said to be TVD if 


ry(„n+i) < ^^(u”). (3.27) 

The following sufficient conditions for (3.22) to be a TVD scheme are due to Harten (1983, 1984): 

TV(R • u”) < rV(«") (3.28a) 


and 


TF(£-u”+‘) >TV(«'‘+*). 


(3.28b) 
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Assuming that the numerical flux in (3.22) is Lipschitz continuous, (3.22) can be written as 

(without the last terms on the implicit and explicit operator, i.e., 2)oa(«j) ) 

- c;_. Vj)”'" = «7+ A(1 - . V. - c;_. A, 

(3.29) 

where are some bounded functions. Then Harten further 

showed that sufficient conditions for (3.29) to be TVD are 

(a) for all i 



(3.30a) 

c+,+c;;.=a(i-9)(c;^.+c;^.)<i, 

(3.30b) 

and 


(b) for all j 


1 

8 

A 

0 
lA 

1 

lA 

o 

(3.31) 

for some finite C. 


To illustrate that (3,22) (using the appropriate numerical flux) satisfies the TVD sufficient 
conditions, we take for example the as Harten’s original numerical flux 


(3.32a) 

with 


Ji = {fi + 9j), 

(3.32b) 

gj = minmod(o-j.+ i Aj._i), 

"i+i = +Ti+i» 

(3.32c) 

(3.32d) 

/ idj+i - 9j)/^j+i Aj.^,1 ^ 0 

^i+j=o, 

(3.32e) 

o-(Sj.+ i) = ^^(Sj.+ i) > 0, 

(3.32f) 

or 

s. . = I ^+5 " 

“'■^1 1 o(“i) ^y+i=0 

(3.32g) 
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Then, for the scheme (3.22) 





Thus, for (3.22) to be TVD, we need 


(3.33) 




= A(1 - 9)C* . 


A(1 - 0) 




+ i>o 


(3.34) 


and 


implies 




(3.35) 




(3.36) 


Therefore, the final scheme using the implicit Euler time discretization (6 = 1) is unconditionally 
TVD and using the trapezoidal method (0 = |) is TVD if lAa^-^ 1 1 < |. 


TVB Schemes : The numerical method (3.22) for an initial-value problem of (3.21) is said to be 
total variation bounded (TVB) in the time interval (0 < < < T) if 


TF(u") < B, (3.37) 

for some fixed S > 0 depending only on u®, all possible n and time step At such that nAt < T. 
TVB schemes are less restrictive than TVD schemes. Clearly TVD implies TVB. There are two 
advantages of TVB schemes over TVD schemes: (a) TVB schemes can be uniformly higher-order 
accurate in space including extrema points; (b) it is easier to devise boundary schemes that are 
TVB for the combined interior and boundary scheme. 


Positive Schemes : The numerical method (3.22) for 0 = 0 for an initial-value problem of (3.21) 
is said to be positive if one can write (3.22) as 

K 

where all the Ck > 0 and = 1. The positive scheme definition for multidimensional 

systems can be found in Liu and Lax (1995). For the proof of whether a higher than first-order 
scheme in multidimensions satisfies the the positivity property, see Einfeldt (1988), and Liu and 
Lax. 
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3.4. Convection Dominating Viscous Flows 


Formal extension of these high-resolution shock-capturing compact schemes to include viscous 
terms while maintaining the same order of accuracy is quite involved and computationally 
expensive (Abarbanel (1994). For convection dominating viscous flows where the accuracy of the 
convection terms is more important, the easiest method for discretizing the viscous terms is the 
standard non-compact second and fourth-order central differencing. Another alternative is to use 
the standard compact fourth-order method. Applying this to yields 



where the C and D operators are defined by 


(3.39a) 


+ Ivi + ^Vi-, , (3.39b) 

(DV)i = V)+, - 2Vj + Vj-, . (3.39c) 

The final scheme using either the non-compact second or fourth-order central differencing or 
(3.39) for the viscous terms is no longer fourth order. How this inconsistent way of discretizing the 
viscous terms affects the overall performance and accuracy of the convection dominating flows 
remain to be addressed. 


IV. Compact High-Resolution Shock-Capturing Schemes 
Based on Total Variation Diminishing in the Mean (TVDM) 

Cockbum & Shu (1994) proposed an explicit compact shock-capturing scheme based on the 
splitting of a Lax-Friedrichs flux and the idea of a total variation diminishing in the mean (TVDM). 
Here, we extend their idea to a family of explicit and implicit schemes with numerical fluxes 
similar to section III but using the TVDM idea. That is, the flux limiters are performed on a 
mean value of U involving a symmetric linear combination of adjacent grid points. However, 
even for scalar hyperbolic conservation laws, their TVDM idea does not completely suppress 
spurious oscillations across discontinuities due to the limiting of the local mean step. See their 
original paper for illustrations. The present straightforward extension of their TVDM idea suffers 
a similar shortcoming. By applying part of their idea and performing flux limiters not on the 
local mean value of U but on U itself, a tremendous improvement in the shock resolution is 
realized. Another alternative in achieving the desired shock resolution is to evaluate the entire 
numerical flux function on the local mean value of U. Before showing that the proposed scheme 
(with the appropriate numerical fluxes) is TVDM for the 1-D scalar case, this high-resolution 
shock-capturing method for the 3-D Euler equations is presented. The reason for presenting the 
schemes in 3-D is to contrast these schemes with the ones proposed in Section III. Although 
the operation count for both families of methods are very close for problems that are lower than 
3-D, for 3-D the TVDM version requires a larger operation count than the Abarbanel and Kumar 
extension. One advantage of the TVDM version is that a spatially and temporally fourth-order 



26 


compact form for both time-accurate and time-marching approaches can be readily obtained as 
opposed to the inability to extend the Abarbanel and Kumar modification in a similar manner for 
time-accurate computations. 

Assume the Gy and in (2.1) are approximated by some compact operator at 
(j Aa;, At/, /Az) 



where for a fourth-order approximation 


(4.1a) 

(4.1b) 

(4.1c) 


i -f- 4J’j-,fc,i -I- , (4.2a) 

^ ^Gj,k+i,i -I- 4Gj,fc,i + Gj,k-i,i ^ , (4.2b) 

{A^H)j^k,i = g , (4.2c) 

^ - Fj-t^k,i ^ , (4.2d) 

{ByG)j,h,i = ^ (^Gj,k+i,i - Gj,k-i,i ^ , (4.2e) 

= I (hj,u,i+i - . (4.2f) 


Similarly, one can define the corresponding 3rd-order and 6th-order compact operators. Since one 
does not gain in operation count for 3rd-order compact schemes over the non-compact cousins 
and 6th-order compact forms are too complicated for higher than 1-D computations (especially for 
3-D), these forms are not discussed here. 

A one-parameter family of explicit and implicit compact schemes suitable for both time-accurate 
and time-marching approaches for (2.1) with 5 = 0 can be written 


= Pm,. -(1 - + M(A-'B,G)2,, + V(A:‘B.B)2J. 


(4.3) 
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Let 


^3,h,i = • (4*4) 

Since A^, Ay, and Az commute, we can multiply the semidiscretized form of (4.3) by A^AyAg 
and (4.3) becomes 


Pill + 

(4.5a) 

A two-parameter family of explicit and implicit counterparts of (4.5a) takes the form 




e 




n-f-1 




i-e 


= + y{A.AyBzH)u 


+ 


03 


1 -|- W 


-f 


(4.5b) 


The one-leg formulation of Dahlquist for this case has the form 


P"tl =P”.»,. - A-(4..1,B.f);, , - \>{A.A.B,G)l, , - X‘(A.A,B,H)':_ , 


+ 


0) 


1 -f- a; 


7t” _ ^ 

^3,h,l ^ j,k,l 


(4.5c) 


The symbol (A) here means the flux evaluations on the right-hand side of (4.5c) are evaluated 
using one of three ways: (1) evaluated atU = f^C^" + (2) evaluated partly at U and 

partly atU = + j^U , and (3) evaluated at U. More details on the various options 

will be discussed shortly. The advantage of (4.5c) over (4.5b), where iterative relaxation methods 
are used to solve the nonlinear algebraic equations as discussed in Sections II and III, carry over 
to the present formulation, especially for the 3-D case. Again, higher than second-order implicit 
counterparts of (4.5b,c) can also be obtained in a similar manner but the resulting scheme involves 
more than three time levels. 

For (4.5) (or (4.4) ) to have a high-resolution shock-capturing property, instead of using (4.2), 
the following is proposed 




(4.6a) 

(4.6b) 

(4.6c) 
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Without lost of generality, it is assumed that U means U for (4.5c) for the following three options 
of evaluating (4.6). In particular, when we say that the flux is evaluated atU = {AgAyA^U)j^k,i, 

we actually mean that the flux is evaluated at {AgAyA^U)j^h,i for (4.5c). 

Option I: To achieve fourth-order spatial differencing, the first option is to evaluate 

j and in exactly the same form and arguments as the and 

in (3.7) - (3.15). This is the author’s first proposed compact form to replace Cockbum 
and Shu’s compact form. The terms in the round brackets on both sides of (4.5) are e.g.. 


(AgAyB^F) ^ ^ - Ag ^j+^,k-i,i) 




= etc... 


(4.7) 


For option I, F* = Fj^r,h+i,i and F* 


= F. 


- 1 1 


, asin n 7'»-n 1SV 


Expanding out the right-hand side of (4.7), one gets 18 terms of the numerical flux evaluations 
in 3-D (in each direction) as opposed to six terms in 2-D (in each direction) and two terms in 
1-D for each of the operators like (4.7) in equation (4.5). Comparing the operation count with 
the Abarbanel and Kumar extension, the TVDM version requires a lot more vector additions 
for the 3-D case. In other words, the vector additions for the Abarbanel and Kumar extension 
increase linearly from 2-D to 3-D but not for the TVDM version. However, the variants of the 
fourth-order Runge-Kutta method of Cockbum and Shu are readily obtained for both time-accurate 
and time-marching approaches as opposed to the inability to extend the Abarbanel and Kumar 
extensions in a similar manner for a spatially and temporeilly fourth-order scheme for time-accurate 
computations. This explicit scheme will be discussed in the next subsection. The viscous analog 
of (4.5) using the fourth-order form (3.39) is straightforward. However, the viscous part is rather 
expensive to compute. The discussion above applies to options II and III as well. 


Option II: The second option is the straightforward extension of Cockbum and Shu’s idea. In this 
case, we define F* i . „ G* , i , and H\ , . i to have the same form as (3.7) - (3.15) except all 

the Rj+i and are evaluated at U U U j+i,h,i and U j+ 2 ,u,i‘ That is, e.g.. 


jp-k 


(4.8) 


where Rj+^ and are the same as in Section III but are evaluated at U. For the corresponding 
MUSCL and flux vector splittings of (3.1), all the limiting is applied to U. 


Preliminary numerical experiments performed by Dr. George Huang of NASA Ames on the 
1-D scalar Burgers equation revealed that the straightforward extension of Cockbum and Shu’s 
idea to the numerical fluxes proposed in section III suffers a shortcoming similar to that of their 
original numerical flux. This TVDM idea does not completely suppress spurious oscillations across 
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discontinuities. However, numerical experiments on option I indicate that spurious oscillations are 
minimized and a tremendous improvement in the shock resolution is realized. 

Option III: The third option is to evaluate all of the terms in the round brackets on the implicit 
and explicit operators in (4.5) at U. That is, e.g., 


E7* 


F + Fj,h,i + 


( 4 . 9 ) 


where Fj+i^h,i, and are the same as in Section III but are evaluated at U. 

For the corresponding MUSCL and flux- vector splitting of (3.1), all the dependent variables and 
limiting are applied to V. Although the third option alters the original compact property of the 
scheme, preliminary numerical experiments show that shock resolution is similar or slightly better 
than option I and far better than option II. In addition, the standard matrix inversion associated 
with the standard compact scheme is not required for option III. In fact, the operation count for 
option III is comparable to the high-resolution case discussed in section III for 2-D, but requires 
more operations for 3-D computations. Note that Option III collapses to the non-compact case for 
1-D unless 5 0 in (2.1a). Although the exact spatial order of accuracy using option III needs 

further investigation, one apparent advantage of option III over their standard non-compact cousins 
is that this scheme can be viewed as an approximation to genuinely multidimensional schemes as 
discussed in Sections II and III. 


Similar to Sections II and III, one can obtain the corresponding ADI delta form for (4.5). For 
steady-state computations, one can use the same simplified first-order spatial differencing for the 
implicit operator (implicit left-hand side) as discussed previously. 

Temporally Higher-Order Explicit TVDM Compact Schemes: Unlike the explicit higher- 

order compact schemes discussed in Section III, the schemes discussed here can retain fourth-order 
spatial and time accuracy for time-accurate and time-marching computations. The explicit scheme 
is the same as (2.15) except the three options in evaluating (4.6) should be applied to the proper 
arguments of each stage of the Runge-Kutta method for high-resolution shock-capturing capability. 
The R° counterpart of (2.15d) is 





( 4 . 10 ) 


Viscous Flows : Formal extension of the TVDM-type compact schemes to include simple viscous 
terms is straightforward using the 4th-order compact operator such as (3.39). However, the grid 
stencil for the inclusion of viscous terms is rather expensive to compute due to the multiplication 
factor of the discretized viscous terms by the Ag^AyAz operator. 

TVDM Property of (^.5) for 1-D Scalar Case: In 1-D with Uj = (Aa«)j, (4.5) becomes 


b7+‘ + A«(B./)7+‘ = 57 - A(1 - »)(B./)7. 


( 4 . 11 ) 
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Cockbum and Shu define the total variation of the mean u by 

oo 

TV(u”)= (4.12) 

J= — OO 

The explicit scheme that Cockbum and Shu considered is for ^ = 0 using the Lax-Friedrichs flux 
splitting. They showed that their explicit scheme satisfies the TVDM sufficient condition, i.e., 

TF(u"+^) < Tr(«”), (4.13) 

under the CFL condition of 

Using the same argument as in Section 3.3 one can readily obtain either the positive or the 
TVDM property for (4.1 1) for the three options of evaluating 


V. Concluding Remarks 

Two families of explicit and implicit compact high-resolution shock-capturing methods for the 
multidimensional compressible Euler equations have been formulated. Some of these schemes can 
be fourth-order accurate away from discontinuities. The attractive property of these compact high- 
resolution shock-capturing schemes is that fourth-order accuracy can be achieved using the same 
grid stencil (5-7 points in each spatial direction) and numerical fluxes as their second-order non- 
compact cousins. In contrast, typical grid stencils for non-compact fourth-order high-resolution 
shock-capturing schemes require 9-11 points in each spatial direction. 

Many variations of these two families of schemes are proposed. The majority of the modified 
Abarbanel and Kumar schemes are best suited for time-marching to the steady state. Although 
they can be used for time-accurate computations, the time accuracy can be at most second order 
in order to have spatially fourth-order accuracy in smooth regions. Two modifications to the 
TVDM idea proposed by Cockbum and Shu to improve shock resolution are discussed. They 
can be used for both time-marching and time-accurate computations. These modifications to the 
scheme of Cockbum and Shu result in far better shock resolution than their original form. In 
2-D, the operation count for both families is comparable, whereas the 3-D TVDM version is more 
expensive to compute at each step than the Abarbanel and Kumar extension. However, a spatially 
and temporally fourth-order compact variant of the Cockbum and Shu scheme is readily obtained 
for both time-accurate and time-marching approaches as opposed to the inability to extend the 
Abarbanel and Kumar modification in a similar manner for time-accurate computations. 

The one-leg formulation of these two families of implicit compact schemes is also proposed. 
If iterative relaxation methods are used, the one-leg forms (for both families) are less expensive 
to compute than their non-one-leg cousins. High-resolution shock-capturing properties of these 
families of compact forms can also be achieved via a variant of the higher-order Lax Friedrichs 
numerical flux. Comparable operations count to their classical shock-capturing counterparts 
can be achieved without the use of Riemann solvers for coupled nonlinear systems. Thus this 
makes high-order high-resolution compact shock-capturing schemes viable and efficient numerical 
methods for 3-D combustion and chemically and thermally nonequilibrium flow computations. 
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Preliminary numerical experiments (performed by Dr. George Huang of NASA Ames) with 
both families of compact schemes for the 1-D scalar Burgers equation and a 2-D shock reflection 
problem indicate the desired shock-capturing property. The time differencing used for the 
numerical experiments for both problems is the implicit Euler method. The Abarbanel and Kumar 
extension gives a better shock resolution than the TVDM version. However, the current TVDM 
version exhibits better shock resolution than the numerical flux used by Cockbum and Shu in their 
original paper. In order to assess the accuracy and shock resolution of these schemes, extensive 
numerical tests for a variety of representative model problems and higher than 1-D Euler equations 
have to be performed. For high gradient, and/or wave propagation and high frequency types of 
flow fields, selection of a limiter with less clipping at extrema points is the key. One possibility 
is the limiter suggested by Davis (1995). The use of limiters for this type of flow containing no 
shock wave can act as a nonlinear method of supplying numerical dissipation. Details of the study 
with limiters other than the ones discussed in section III are the subject of current research. 

Without extensive numerical experiments in hand, all we can comment on is the relative 
advantages and disadvantages of these two families of schemes based on operations count (vector 
additions and tridiagonal matrix inversion due to the compact formula), the requirement of a 
special numerical boundary treatment (Carpenter et al. 1993), applicability to and accuracy 
for time-accurate calculations, and ease of implementation to include source terms. From the 
discussion in previous sections, one can conclude that for all but the fourth-order time-accurate 
capability issue, the Abarbanel and Kumar extension appears to be more efficient than the TVDM 
version. 

The majority of these implicit schemes are especially suited for time-marching approaches 
to steady-state numerical solutions, since higher-order spatial accuracy can be achieved with 
minimal effort and the steady states are independent of the time step. The proper choice of time 
discretizations for the proposed two families of schemes for wave propagation, and computational 
aeroacoustics type of applications has not been addressed. Thus, the time discretizations discussed 
in this paper might not be optimized for the particular types of flows in the sense of wave 
resolution and phase error discussed in Tam and Webb, Lele, Davis, and the 1995 Workshop on 
Aeroacoustics. This topic is ongoing research. Future study will include practical application of 
these families of methods to a variety of flow physics. 
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